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CONSISTENCY  TESTS  OF  ACOUSTIC  PROPAGATION  MODELS 

by 

Finn  B.  Jensen  and  William  A.  Kuperman 


I  ABSTRACT 

Three  wave-theory  models  (NM,  FFP,  PE)  and  one  ray  model  are  applied  to 
four  different  ocean  environments:  a  range-dependent  surface  duct,  a 
deep-water  environment  with  a  homogeneous  bottom,  a  shallow-water 
environment  with  a  homogeneous  bottom,  and  a  si  oping- bottom  environment 
with  a  layered  bottom.  The  consistency  among  the  acoustic  models  is 
clearly  demonstrated  through  the  agreement  between  model  results  for  the 
various  test  problems. 
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INTRODUCTION 


In  this  report  we  present  the  results  of  an  inter-model  comparison  for 
which  four  different  underwater  acoustic  environments  were  studied.  Three 
wave-theory  models  were  used:  Normal  Mode  (NM),  Fast  Field  Program  (FFP), 
and  Parabolic  Equation  (PE).  This  study  is  an  outgrowth  of  the  Parabolic 
Equation  Workshop  <1>  held  at  the  US  Naval  Ocean  Research  and  Development 
Activity  (NORDA)  in  April  1981.  The  test  cases  were  provided  by  the 
Numerical  Modelling  Division  of  NORDA  and  various  participants  used  these 
test  cases  for  their  PE  models;  comparisons  between  PE  results  were 
presented  at  the  Workshop.  We  also  did  the  test  cases  with  the  PE,  but  in 
addition  we  ran  most  of  the  cases  with  other  wave  models  in  order  to  obtain 
independent  solutions  to  the  four  problems.  On  occasion  we  even  ran  a  ray 
model  for  completeness.  The  result  is  a  rather  complete  set  of  test  cases, 
each  of  which  has  been  studied  with  various  models.  Where  there  is 
essentially  full  agreement  among  the  three  wave  models  for  the  complicated 
environments  (which  do  not  lend  themselves  to  analytic  solutions),  we 
believe  that  the  answers  presented  are  the  correct  solutions  to  the  wave 
equation  since  the  numerical  algorithms  used  are  totally  different  for  the 
three  models  (and  the  PE  model  even  uses  an  approximation  of  the  wave 
equation! ). 

In  the  next  four  chapters  we  present  the  NORDA  test  cases  together  with 
numerical  results  from  the  various  models.  We  then  summarize  the  results 
of  this  inter-model  comparison  and  give  some  numerical  parameters  and 
associated  computer  times  for  NM  and  PE  solutions  to  the  four  test 
problems.  Appendix  A  outlines  the  mathematical  basis  of  the  various 
acoustic  models  and  Appendix  B  provides  a  brief  description  (with 
appropriate  references)  of  the  specific  computer  codes  used. 
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1  TEST  CASE  1:  RANGE-DEPENDENT  SURFACE  DUCT 

This  test  problem  has  a  range-dependent  sound  speed  and  consists  of  profile 
1A  from  0  to  20  km,  a  transition  from  profile  1A  to  IB  over  the  range 
interval  20  to  30  km,  and  profile  IB  from  30  to  50  km,  as  shown  in  Table  1. 

Figures  1,  2,  and  3*  display  the  result  for  the  three  receiver  depths.  The 
NM  and  FFP  solutions  are  applicable  only  in  the  range- independent  region, 
that  is,  up  to  the  dashed  line  at  range  20  km.  Beyond  20  km  the  difference 
between  the  FFP  and  PE  solutions  indicates  the  effect  of  range  dependence 
on  computed  propagation  losses.  We  see  that  at  ranges  of  less  than  20  km 
there  is  excellent  agreement  between  the  results  of  the  three  models,  the 
greatest  disagreement  occurring  in  the  nearfield  where  the  FFP  should  work 
best  (see  App.  A).  Note  though,  that  in  the  nearfield  the  PE  and  FFP  have 
better  agreement  with  each  other  than  they  do  with  the  NM  calculation.  Of 
course  the  NM  calculation  does  not  include  the  nearfield  contribution,  but 
it  is  interesting  to  note  how  well  the  PE  performs  in  the  nearfield. 

The  main  difficulty  in  running  this  problem  with  the  PE  model  was  to  avoid 
acoustic  returns  from  the  lower  boundary  of  the  truncated  z-domain.  We 
obtained  stable  results  by  moving  the  boundary  to  a  depth  of  4000  m,  and  by 
introducing  an  artificial  absorption  in  the  lower  1000  m  of  that  domain. 
To  obtain  discrete  normal-mode  solutions  for  this  case,  which  essentially 
is  a  "continuous  spectrum"  problem,  we  had  to  introduce  an  artificial 
bottom  with  a  speed  of  1500  m/s  (equal  to  maximum  speed  in  the  surface 
duct).  By  moving  this  bottom  to  a  depth  of  4000  m  we  obtained  stable  NM 
solutions  to  cases  1A  and  IB,  with  the  trapped  "continuous  mode"  in  the 
surface  duct  now  being  simulated  through  a  set  of  36  discrete  modes. 
However,  for  the  receiver  below  the  surface  channel  (case  1C)  we  found  no 
efficient  manner  in  which  the  discrete  NM  model  could  simulate  the 
"continuous  spectrum"  problem. 

The  conclusion  from  Figs.  1  to  3  is  that  the  PE  gives  the  correct  result  in 
the  range-independent  part  of  the  environment  and  because  of  the 
mathematical  nature  of  the  PE  solution  (see  App.  A)  there  is  no  reason  to 
believe  that  the  accuracy  of  this  solution  diminishes  in  the  mildly  range- 
dependent  region. 

Contouring  the  PE  transmission-loss  results  over  depth  and  range  provides 
further  insight  into  the  basic  propagation  mechanisms  involved  in  this 
case.  Figure  4  shows  propagation  from  the  deep  to  the  shallow  surface  duct 
with  the  source  in  the  duct  (25  m).  The  contour  lines  indicate  that  there 
is  essentially  one  mode  propagating  in  the  duct  out  to  ranges  of  about 
25  km.  Here  the  mode  is  cut  off  and  sound  is  leaving  the  duct  as  indicated 
by  the  dashed  arrow.  Such  mode  cut-off  and  radiation  into  the  bottom  has 
been  studied  in  detail  elsewhere  <2>.  If  we  reverse  the  propagation 
direction  we  get  the  result  shown  in  Fig.  5.  Here  we  have  placed  the 
source  below  the  shallow  duct  (SD  =  250  m).  We  see  that  sound  is 
propagating  with  high  losses  in  the  duct  out  to  ranges  of  about  25  km. 
Here  the  duct  widens  and  sound  becomes  trapped  in  the  duct;  consequently 
propagation  conditions  improve  (longer  spacing  between  contour  lines). 


All  figures  are  at  the  end  of  the  text. 
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TABLE  1 


TEST  CASE  1 


Frequency  =  25  Hz 

Source  depth  (SO)  =  25  m 

Receiver  depths  (RD)  =  25,  250  and  400  m 


Max 

Depth 

range  = 

Sound  speed 
(m/s) 

1480 

50  km 

Profile  1A 

Density 

Attenuation 

(m) 

0 

(g/cnr*) 

1. 

(dB/km) 

0. 

300 

1500 

1. 

0. 

1000 

1460 

1. 

0. 

OO 

1460 

1. 

0. 

Depth 

Sound  speed 

Profile  IB 

Density 

Attenuation 

(m) 

(m/s) 

(g/cm3) 

(dB/km) 

0 

1492.38 

1. 

0. 

200 

1505.71 

1. 

0. 

300 

1500.00 

1. 

0. 

1000 

1460.00 

1. 

0. 

OO 

1460.00 

1. 

0. 

The  profile  in  the  transition  interval  is  range  dependent  between  depths  of 
0  and  300  m  and  is  given  by: 

h(r)  =  300  -  10  (r-20) 

Cs(r)  =  1480  +  1.238  (r-20) 

Cd(r)  =  1500  +  0.571  (r-20), 

where 

h(r)  is  the  depth  of  the  duct  in  metres 
Cg(r)  is  the  sound  speed  (in  m/s)  at  zero  depth 
Cd(r)  is  the  sound  speed  (in  m/s)  at  the  bottom  of  the  duct 
r  is  horizontal  range  in  km. 

Plots: 

Separate  plots  from  0  to  50  km  for  each  of  the  three  receiver  depths. 
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ill 


The  results  in  Figs.  4  and  5  naturally  lead  to  a  test  for  reciprocity  in 
the  PE  approximation  to  the  wave  equation.  The  result  for  case  IB  is  shown 
in  Fig.  6;  the  reciprocal  case  for  the  deep  to  shallow  with  SD  =  25  m  and 
RD  =250  m  is  shallow  to  deep  and  SD  =  250  m  and  RD  =  25  m.  That  is,  not 
only  is  the  direction  of  propagation  reversed  but  the  source  and  receiver 
depths  must  also  be  interchanged.  Reciprocity  for  range-dependent 
propagation  requires  that  only  at  the  final  range  of  propagation  must  the 
answers  be  the  same,  and  indeed  they  are  at  a  range  of  50  km.  The 
calculations  are  not  the  same  throughout  the  total  range  for  a 
range-dependent  environment  because  the  environments  are  different  for  the 
two  cases.  For  example,  at  20  km,  the  answers  are  far  apart  because  in  one 
case  we  are  propagating  in  the  shallow  duct  (no  trapping  of  sound)  and  in 
the  other  case,  in  the  deep  duct  (one  mode  trapped).  Figure  6,  however, 
demonstrates  that  reciprocity  is  included  in  the  parabolic-equation  method. 


2  TEST  CASE  2: 

RANGE- INDEPENDENT  ENVIRONMENT  WITH  DIFFERENT  BOTTOM  SPEEDS 

This  test  problem  has  a  range- independent  environment  with  a  bi -linear 
sound  speed  in  the  water  column  and  a  half-space  bottom.  It  consists  of 
three  parts,  A,  B  and  C,  with  different  sound-speed  profiles,  as  shown  in 
Table  2. 

TABLE  2 
TEST  CASE  2 
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Since  this  case  is  a  range- independent  environment  we  know  that  NM  and  FFP 
are  applicable,  with  the  caveat  that  the  discrete  NM  solution  will  not  give 
an  adequate  description  of  the  nearfield  (see  App.  A).  The  results  for 
the  three  cases  are  shown  in  Figs.  7  to  9.  Note  first  the  excellent 
agreement  between  FFP  and  NM,  indicating  that  we  have  the  true  solution  to 
the  wave  equation.  Next  notice  that  the  agreement  between  PE  and  NM/FFP 
deteriorates  with  increasing  bottom  speed  even  though  the  approximate 
levels  are  correct.  With  increasing  bottom  speed  the  maximum  angle  of 
propagation  (with  respect  to  the  horizontal)  increases  and  so  does  the 
number  of  propagating  modes,  and  both  factors  affect  the  accuracy  of  the  PE 
solution.  For  cases  A,  B,  and  C  the  angles  are  15°,  30°,  and  40° 
respectively,  and  the  number  of  modes  are  11,  22,  and  28. 

The  PE  is  a  narrow-angle  approximation  to  the  wave  equation  (see  App.  A) 
and  accurate  results  are  obtained  only  for  propagation  angles  of  less  than 
15  to  20°  with  the  horizontal.  Moreover,  there  is  an  inherent  phase  error 
in  the  PE  that  results  in  slightly  different  spatial  interference  patterns 
than  those  calculated  by  an  exact  solution  to  the  wave  equation  (NM  or 
FFP).  This  phase  error  becomes  more  evident  as  the  number  of  modes 
increases.  In  this  study  we  have  employed  a  correction  procedure  <3>  that 
should  reduce  PE  phase  errors,  particularly  in  deep-water  type  environ¬ 
ments.  For  case  2A  (Fig.  7),  where  propagation  angles  are  less  than  15°, 
the  PE  follows  the  correct  interference  pattern  quite  accurately, 
indicating  that  the  phase-correction  procedure  works  properly.  However, 
with  increasing  angle  of  propagation  (Figs.  8  and  9)  the  accuracy  of  the  PE 
solution  deteriorates,  and  though  some  prominent  features  are  retained 
(including  a  correct  mean  level),  a  component  of  randomness  is  seen  in  the 
PE  interference  pattern  for  these  wide-angle  cases. 

We  showed  for  case  1  (Ch.  1)  that  the  PE  includes  the  "continuous  spectrum" 
contribution.  In  Appendix  A  we  discuss  this  aspect  in  more  detail  with 
reference  to  case  2A,  showing  that  the  PE  substantially  agrees  with  the 
nearfield  FFP  result,  the  classical  Lloyd  mirror  pattern. 


3  TEST  CASE  3: 

RANGE- INDEPENDENT  SHALLOW-WATER  ENVIRONMENT 


This  test  problem  has  a  range- independent  environment  and  consists  of  an 
isovelocity  water  column  over  an  isovelocity  half-space  bottom.  There  are 
two  parts,  A  and  B,  with  different  source  and  receiver  depths,  as  shown  in 
Table  3. 

Again  we  see  that  for  both  of  these  cases  (Figs.  10  and  11)  we  have 
agreement  between  FFP  and  NM  and  hence  we  have  the  wave-theory  solution. 
For  this  environment  the  steepest  angle  of  propagation  is  around  20°  and 
hence  the  narrow-angle  approximation  should  not  affect  the  accuracy  of  the 
PE  solution,  at  least  not  for  source  and  receiver  at  mid-depth.  Since  the 
PE  phase-correction  procedure  <3>  is  not  effective  for  this  environment 
with  heavy  bottom  interaction,  we  have  tried  to  minimize  phase  errors  in 
the  PE  solution  by  an  appropriate  choice  of  reference  speed  (see  App.  A). 
The  reference  speed  should  be  chosen  to  match  the  phase  velocity  of  the 
most  important  mode  <4>.  For  case  3A,  with  source  and  receiver  at 
mid-depth,  sound  is  mainly  propagating  in  the  lower-order  modes  (the  total 
number  of  modes  is  11).  By  choosing  a  PE  reference  speed  of  1500  m/s, 
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which  is  close  to  the  phase  velocity  of  the  first  mode  (1500.6  m/s),  we 
obtained  the  good  agreement  between  PE  and  NM/FFP  solutions  shown  in 
Fig.  10. 

For  case  3B,  with  source  and  receiver  near  the  bottom  boundary,  higher- 
order  modes  are  significantly  excited,  making  the  interference  pattern  more 
complex  than  for  case  3A.  We  see  in  Fig.  11  that  the  PE  solution  does  not 
precisely  track  the  true  interference  pattern.  Nevertheless,  observe  the 
agreement  at  the  large  destructive  interference  null  at  a  range  of  7  km. 
This  agreement  was  obtained  using  a  PE  reference  speed  of  1550  m/s,  which 
is  a  mean  value  for  the  phase  speeds  of  this  problem.  (The  difference 
between  FFP  and  NM  solutions  around  7  km  is  attributed  to  inaccuracies  in 
the  computations;  phase  accuracy  is  crucial  for  achieving  the  degree  of 
cancellation  between  modes  shown  here).  Note  in  Fig.  11  that  the  average 
level  for  the  PE  solution  is  a  little  too  low.  This  is  probably  due  to  the 
narrow-angle  approximation,  which  affects  propagation  amplitudes  at  steeper 
angles. 

The  most  interesting  feature  of  this  test  problem  is  the  large  destructive 
interference  null  for  case  3B.  The  transmission-loss  contours  of  Fig.  12 
created  by  the  PE  model  show  the  spatial  extent  of  the  null  region 
around  7  km.  A  special  check  on  the  adequacy  of  the  gauss i an  beam  for 
initializing  the  PE  solution  was  done  for  case  3;  we  obtained  identical 
results  whether  starting  with  a  normal-mode  field  or  with  the  gaussian 
beam,  even  for  the  source  close  to  the  bottom. 


TABLE  3 
TEST  CASE  3 


Frequency 

— 

250  Hz 

Water  depth 

— 

100  m 

Sound  speed  in  water 

“ 

1500  m/s 

Density  in  water 

= 

1.0  g/cm3 

Density  in  bottom 

- 

1.2  g/cm3 

Attenuation  in  water 

- 

0 

Bottom  attenuation 

— 

0.5  dB/\ 

Bottom  sound  speed 

= 

1590  m/s 

Max  range 

10  km 

Part  A: 

Part  B 

Source  depth 

50  m 

99.5  m 

Receiver  depth 

50  m 

99.5  m 

Plots: 

Separate  plots  from  5  to 

10 

km  for  each  of  the 

two  parts. 
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4  TEST  CASE  4: 

PROPAGATION  OVER  A  SLOPING  BOTTOM  WITH  A  SEDIMENT  LAYER 


This  problem  has  a  range- independent  sound  speed  in  the  water  column,  a 
range-dependent  bottom  depth  and  a  simple  geoacoustic  bottom  that  is  tied 
to  the  sound  speed  at  the  bottom  of  the  water  column,  as  shown  in  Table  4. 

This  is  a  fairly  complicated  environment,  where  the  steepest  propagation 
angles  are  around  40°.  Hence  the  narrow-angle  approximation  could  cause 
inaccuracies  in  the  PE  results,  particularly  at  short  ranges  where  steep 
bottom-bounce  paths  are  most  important.  No  phase  correction  was  applied  to 
the  PE  solution;  we  chose  a  reference  speed  of  1500  m/s.  The  total  number 
of  modes  for  this  case  is  77.  The  range-averaged  propagation-loss  curves 
were  produced  using  a  sliding  rectangular  window  (equal  weighting  of  all 
samples). 

Figure  13  shows  the  bathymetry,  the  layering,  and  a  ray  trace  of  the 
waterborne  paths  to  indicate  the  type  of  long-range  propagation  expected. 
Figures  14  and  15  show  good  agreement  in  the  levels  obtained  by  the  NM  and 
FFP  calculations.  In  particular,  the  20  km  smoothed  results  (Fig.  15)  are 
virtually  identical  for  FFP  and  NM.  Figures  14  and  15  also  indicate  good 
agreement  with  the  PE  solution.  The  slightly  lower  PE  levels  for  ranges 
less  than  25  km  are  probably  caused  by  an  incorrect  handling  of  bottom- 
bounce  paths  (not  traced  in  Fig.  13  to  keep  the  plot  simple).  These  steep 
propagation  paths  will  carry  only  little  energy  after  a  few  bounces  off  the 
bottom.  Note  in  Fig.  15  that  NM  and  FFP  solutions  are  applicable  only  out 
to  a  range  of  150  km.  Beyond  this  range  we  have  bathymetric  changes  that 
are  included  only  in  the  PE  solution.  We  see  that  for  this  particular 
receiver  depth  there  is  no  effect  of  the  changing  water  depth  on 
propagation  loss  between  150  and  190  km.  Then  there  is  a  sharp  fall -off  in 
the  PE  curve  as  the  receiver  moves  into  the  bottom. 

Figure  16  compares  the  PE  result  with  range-dependent  ray  theory  and 
adiabatic  mode  theory.  Neither  the  adiabatic  mode  nor  the  ray  calculation 
could  be  made  for  the  full  range,  since  these  particular  models  were  not 
set  up  to  have  a  receiver  in  the  bottom.  For  the  ray-theory  calculation 
the  sediment  was  taken  to  have  a  density  of  one  and  no  attenuation.  This 
was  done  so  that  the  sediment  layer  could  be  treated  as  part  of  the  sound- 
speed  profile  in  the  water  column.  In  general  there  is  good  agreement 
among  the  models.  The  slightly  higher  level  of  the  ray-theory  result 
beyond  75  km  is  because  no  attenuation  was  included  in  the  sediment. 

Results  for  test  case  4C  are  shown  in  Fig.  17,  where  PE  is  compared  against 
adiabatic  mode  theory.  There  is  agreement  between  levels,  but,  as  we 
expected,  the  interference  pattern  is  slightly  shifted  because  of  phase 
errors.  Figure  18  displays  test  4D  results  and  compares  NM,  FFP  and  PE. 
Note  that  ranges  less  than  150  km  correspond  to  a  range- independent 
environment  where  NM  and  FFP  are  applicable.  Beyond  150  km  comparison  of 
the  PE  result  with  the  NM/FFP  results  indicates  the  difference  between 
constant  depth  and  sloping-bottom  propagation.  The  initial  increase  in 
level  up  the  slope  is  a  geometrical  effect  of  the  decrease  in  water  depth. 
In  the  shallow-water  region  beyond  200  km  the  increased  interaction  with 
the  bottom  results  in  increased  losses  with  range. 

Finally,  Fig.  19  indicates  that  the  PE  result  for  this  case  is  in  good 
agreement  with  adiabatic  mode  theory  while  ray  theory  again  gives  levels 
that  are  too  high  because  no  attenuation  was  included  in  the  sediment  layer. 
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TABLE  4 

TEST  CASE  4 


Frequency 
Source  depth 
Receiver  depths 
Max  range 

Water  column  parameters 


=  2S  Hz 
=  600  m 

=  150  and  700  m 
=  250  km 


Depth 

Sound  speed 
(m/s) 

Density 

Attenuation 

(m) 

(g/cm3) 

(dB/km) 

0 

1539.3 

1. 

0. 

30 

1539.8 

1. 

0. 

200 

1534.2 

1. 

0. 

600 

1502.4 

1. 

0. 

700 

1495.4 

1. 

0. 

800 

1491.8 

1. 

0. 

1000 

1488.9 

1. 

0. 

1100 

1487.5 

1. 

0. 

1200 

1487.9 

1. 

0. 

3410 

1525.0 

1. 

0. 

Bottom  depth 

Ranqe 

Bottom  depth 

0-150  km 

Constant  3410  m 

150-200  km 

linearly  decreasing  from 

3410  m  to  200  m  with 

increasing  range 

200-250  km 

Constant  200  m 

Sediment: 

Sound  speed  at  top  of  sediment  is: 

CTs(r)  =  0.975  x  Cw(r), 

where  C  ( r)  is  the  sound  speed  at  the  bottom  of  the  water  column. 

Sound  speed  at  the  bottom  of  sediment  is: 

CB$(r)  =  1.305  x  CTs(r) 

Sediment  thickness  =  454  m  (constant  in  range) 

Attenuation  =  0.375  dB/km  =  0.015  dB/km/Hz 
Density  =1.5  g/cm3 

Basement: 

Sound  speed  =  C85(r) 

Attenuation  =  0. 

Density  =2.5  g/cm3 

Example:  At  a  range  of  200  km, 

Water  depth  =  200  m 

Cw(200  km)  =  1534.2  m/s 
CT^(200  km)  =  1495.8  m/s 
CBS(200  km)  =  1952. 1  m/s 
Basement  depth  =  654  m 
Plots: 

For  700  m  receiver  depth: 

case  A  -  1  km  average  from  10  to  50  km 
case  B  -  20  km  average  from  0  to  250  km 

For  150  m  receiver  depth: 

case  C  -  1  km  average  from  175  to  225  km 
case  0  -  20  km  average  from  0  to  250  km 
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SUMMARY  AND  CONCLUSIONS 

Considering  the  completeness  of  this  study  (four  models  compared)  and  the 
realistic  ocean  environmental  descriptions  contained  in  the  four  test 
problems,  it  seems  likely  that  the  results  given  in  this  report  could 
constitute  a  set  of  benchmark  tests  for  other  propagation  models  in  the 
future.  With  that  prospect  in  view  we  have  attempted  to  include  as  much 
useful  information  as  possible  on  the  various  test  cases  and  on  the 
acoustic  models.  In  this  spirit  and  to  facilitate  the  reproduction  of  the 
results  in  this  report,  we  have  listed  in  Table  5  representative  numerical 
parameters  and  associated  computer  times  for  both  the  normal  mode  and  the 
PE  model.  The  necessary  sampling  in  depth,  Az,  for  obtaining  stable 
numerical  solutions  is  seen  to  be  of  the  same  order  of  magnitude  for  both 
models.  The  sampling  in  range,  Ar,  has  significance  for  numerical  accuracy 
only  in  the  PE  model.  The  transform  size  given  for  each  test  problem  can 
be  used  to  compute  the  total  „ depth  domain  (H  )  used  in  the  PE 
computation,  i.e.  Hmax  =  Az  •  2  .  The  computer3  times  relate  to  a 

UNIVAC  1100/60.  We  see  that  the  PE  model  is  particularly  slow  when  there 
is  heavy  bottom  interaction  (cases  3  and  4).  We  did  not  attempt  to 
optimize  the  FFP  runs.  However,  an  FFP  calculation  is  generally  slower 
than  a  mode  computation. 

TABLE  5 

NUMERICAL  PARAMETERS  ANO  ASSOCIATED  COMPUTER  TIMES  (UNIVAC  1100/60) 


FOR  MODE  AND  PE  SOLUTIONS  TO  THE  FOUR  TEST  CASES 
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To  sum  up,  we  have  presented  numerical  calculations  for  the  NORDA  test 
cases.  The  consistency  among  normal  mode,  FFP,  and  PE  (and  ray  theory)  has 
been  clearly  demonstrated  and  cases  of  disagreement  are  well  understood.  It 
was  never  a  straightforward  procedure  to  obtain  correct  answers  from  either 
of  the  models  employed  in  this  study.  Unfortunately  the  state-of- 
development  of  scientific  propagation  models  is  such  that  considerable 
experience  with  a  given  model  is  necessary  in  order  to  get  meaningful 
results.  However,  of  the  models  used  here  (see  App.  B),  the  ray  and 
normal-mode  programs  definitely  have  reached  a  higher  level  of  refinement 
and  automation  than  the  other  two  models  and  they  are  therefore  the  easiest 
to  run.  The  PE  model  and  particularly  the  FFP  model  require  at  present 
considerable  insight  into  the  specific  numerical  schemes  in  order  to 
produce  meaningful  propagation-loss  predictions. 
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FIG .  1  INTER-MODEL  COMPARISON  FOR  CASE  1A 
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FIG.  2  INTER-MODEL  COMPARISON  FOR  CASE  IB 
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FIG.  3  INTER-MODEL  COMPARISON  FOR  CASE  1C 
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FIG.  4  CONTOURED  PROPAGATION  LOSS  FROM  PE  MODEL  GOING  FROM  DEEP  TO  SHALLOW  DUCT 
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SACLANTCEN.  TEST  IB  —  SHALLOW/DEEP 


FIG.  5  CONTOURED  PROPAGATION  LOSS  FROM  PE  MODEL  GOING  FROM  SHALLOW  TO  DEEP  DUCT 
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FIG.  6  RECIPROCITY  TEST  FOR  PE  MODEL  FOR  CASE  IB 
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SACLANTCEN,  TEST  2A 

FIG.  7  INTER-MODEL  COMPARISON  FOR  CASE  2 A 
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FIG .  8  INTER-MODEL  COMPARISON  FOR  CASE  2B 
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FIG.  9  INTER-MODEL  COMPARISON  FOR  CASE  2C 
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FIG.  10  INTER-MODEL  COMPARISON  FOR  CASE  3 A 
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FIG .  11  INTER-MODEL  COMPARISON  FOR  CASE  3B 
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SACLANTCEN.  TEST  3B 

FIG.  12  CONTOURED  PROPAGATION  LOSS  FROM  PE  MODEL  FOR  CASE  3B 
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FIG.  13  RAY  TRACE  FOR  CASE  4.  The  dashed  line  indicates  the  water/sediment  interface 
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FIG.  14  INTER-MODEL  COMPARISON  FOR  CASE  4A 
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FIG.  15  INTER-MODEL  COMPARISON  FOR  CASE  4B 
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FIG  ,  16  COMPARISON  OF  THREE  RANGE-DEPENDENT  MODEL  RESULTS  FOR  CASE  4B 
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FIG.  17  INTER-MODEL  COMPARISON  FOR  CASE  4C 
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FIG.  18  INTER-MODEL  COMPARISON  FOR  CASE  4D 
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FIG.  19  COMPARISON  OF  THREE  RANGE-DEPENDENT  MODEL  RESULTS  FOR  CASE  4D 
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APPENDIX  A 


MATHEMATICAL  FOUNDATION  OF  ACOUSTIC  MODELS 


In  this  appendix  we  briefly  present  the  mathematical  foundations  of  the 
models  used  in  the  main  body  of  this  report.  A  more  detailed  description 
can  be  found  in  references  A.l  to  A. 9.  Appendix  B  describes  the  specific 
models  used. 

The  starting  point  for  all  the  models  is  the  wave  equation  for  a  harmonic 
point  source  with  time  dependence  exp(-iuit), 


V2«|»(x,y,z)  + 


r  i 

c(x,y,z)  J 


<Kx,y,z)  =  -6(x-xo)6(y-yo)6(z-zo) 


»|)  =  <]>exp(-iu)t) 


(Eq.  A.l) 
(Eq.  A. 2) 


At  any  point  (x,y,z)  in  the  medium,  the  velocity  potential  $  satisfies 
Eq.  A.l,  where  c(x,y,z)  is  the  sound  speed  of  the  medium  and  6  is  the  Dirac 
delta  function.  The  source  is  at  the  coordinate  (x  ,y  ,z  ),  where  z 
is  the  depth  coordinate,  which  is  taken  to  be  positivi  ui  the  downward 
direction  from  the  ocean  surface. 

For  the  boundary  condition  at  the  ocean  surface  we  take  the  density  of  air 
to  be  negligible  compared  with  that  of  water;  hence,  the  pressure  must 
vanish  at  the  ocean  surface  ("pressure-release  surface").  At  a  boundary 
between  two  media  such  as  the  ocean  and  the  ocean  bottom,  the  balancing  of 
forces  at  the  interface  require  that  physical  quantities  such  as  particle 
velocity  and  pressure  be  continuous  across  the  boundary: 


_  9<|>  . 

"  3xi  * 


p  =  -  iwp<J) 


xi  =  x,y,  or  z 


(Eq.  A. 3) 


(Eq.  A. 4) 


If  the  ocean  bottom  is  treated  as  an  elastic  medium  that  can  support  shear 
motions,  there  is  the  additional  boundary  condition  that  tangential  stress 
must  be  continuous.  Since  the  water  column  cannot  support  shear  waves, 
this  requires  that  the  tangential  stress  in  the  ocean  bottom  vanishes  at 
the  interface. 

The  four  widely  used  solution  techniques  for  Eq.  A.l  are  schematically 
represented  below  in  Fig.  A.l. 


HAVE  EQUATION 


high-frequency 

approximation 


FAST  FIELD 


range*- independent 
wave  theory 


NORMAL  MODE 


range-dependent 
wave  theory 


PARABOLIC- EQ. 


FIG.  A.l  TECHNIQUES  FOR  SOLVING  THE  WAVE  EQUATION 
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We  now  briefly  describe  the  derivation  of  these  solutions,  starting  with 
range- independent  wave  theory. 


A.  1  Fast  Field  Solution 

Here  we  are  solving  the  wave  equation  for  the  case  where  the  sound-speed 
profile  is  only  a  function  of  depth  and  the  bottom  is  flat;  this  type  of 
environment  is  often  referred  to  as  the  horizontally  stratified  ocean.  From 
Eq.  A.  1  we  therefore  have  that  c(x,y,z)  is  simply  c(z).  Because  the 
environment  is  independent  of  "r",  the  horizontal  coordinates  (x,y),  one 
possible  method  of  solving  Eq.  A.l  is  to  Fourier  decompose  the  acoustic 
field  into  an  infinite  set  of  horizontal  waves: 

<l>(x,y,z)  =  — /  d2q  u(q,z)  ein‘r  .  (Eq.  A. 5) 


Substituting  Eq.  A.  5  into  Eq.  A.l  we  obtain  the  equation  for  <Knx»ny  ;z) , 
a2u(nx,nv;z)  , 

- ^£2* —  +  [k2(z)  -  n2Mnx.ny;z)  =  -  ^  &(z  -  zQ),  (Eq.  a. 6) 


where  k(z)  =  u»/c(z)  and  n2  =  q*  +  the  horizontal  wavenumber  of  the 

x  y 

individual  plane  waves. 

Using  polar  coordinates  we  can  rewrite  Eq.  A. 5  as 

♦(r.z)  =h  j  d0  j  n  dn  u(n,z)e1nrcos6  .  (Eq.  A. 7) 


We  now  integrate  over  the  azimuthal  angle  to  obtain 

<Kr,z)  =  j  n  dq  u(n,z)JQ(nr)  ,  (Eq.  A. 8) 

where  JQ  is  the  zeroth  order  Bessel  function.  Using  the  relationship  that 
J0(nr)  =  \  tHo1)(nO  +  H^2)(nr)]  , 


where  the  H's  are  Hankel  functions  and  noting  from  Eq.  A. 6  that  u(n,z)  is 
even  in  q»  we  can  rewrite  Eq.  A. 8  as 


<J>(r 


= ?  L 


n  dq  u(n,z)  , 


(Eq.  A. 9) 
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m 


where  now  the  integration  over  q  is  from  -»  to  ».  For  ranges  greater  than 
a  few  wavelengths  from  the  source,  the  asymptotic  form  of  the  Hankel 
function  can  be  used: 


H<1}(qr)  ~  (2/nr)r)'s  exp[i(qr  -  n/A )] 
and,  hence,  Eg.  A. 9  can  be  expressed  as 


4>(r,z)  = 


1_  C° 

42.  n  4rJ-° 


dn  Vnu(n,z)eir|r  , 


(Eq.  A. 10) 


where  the  factor  1/Vr  indicates  cylindrical  spreading. 


Equation  A. 10  can  be  numerically  integrated  to  obtain  the  acoustic  field  at 
the  range  r  and  depth  z.  In  order  to  do  this  we  must  solve  Eq.  A. 6  for 
many  q's  to  have  a  sufficient  set  of  u‘s  as  a  function  of  q  so  that  the 
integration  over  q  in  Eq.  A. 10  can  be  performed.  Given  that  u  has  has  been 
obtained  numerically  as  a  function  of  q,  the  integration  can  be  done  using 
an  FFT  algorithm.  This  total  procedure  is  called  the  Fast  Field  Program 
(FFP)  <A.10,  A.ll>,  although  most  of  the  numerical  effort  goes  into  solving 
Eq.  A. 6  for  the  many  q‘s.  For  computation  we  discretize  Eq.  A. 10  by 
letting 


0=0+  mAq ;  r  =  r  +  nAr;  (m,n)  =  0,1,2...,N  -  1  , 


with  the  additional  relation 


(Eq.  A. 11) 


ArAq  = 


(Eq.  A. 12) 


and  N  is  an  integral  power  of  two.  Note  that  the  discretization  relations 
of  Eq.  A.  11  restricts  the  solution  to  outgoing  waves.  Substituting 
Eq.  A. 11  into  Eq.  A. 10  we  obtain 


-in/4 


4>(r  z)  =  Aq  ■  ■ 
n  V2nr 


iq  r 
e  '0  n 


N-l 

£  Xme 

M=0  m 


i2nmn/N 


(Eq.  A. 13) 


and  hence  the  input  to  the  FFT  is 


X  =  Vnm  u(0m,z)eimVn 
m  'm  'm 


(Eq.  A. 14) 


Equations  A.  13  and  A.  14  specify  the  numerical  procedure  to  be  employed  in 
solving  the  wave  equation  using  the  FFP  approach  after  Eq.  A. 6  has  been 
solved  numerically  for  the  complete  set  of  u's  as  a  function  of  q.  As 
mentioned  above,  the  main  effort  in  the  FFP  approach  is  the  numerical 
integration  of  Eq.  A. 6  and  not  the  final  implementation  of  Eqs.  A.  13  and 
A.  14,  which  is  a  simple  FFT  computation.  Numerical  procedures  for 
integrating  Eq.  A. 6  are  given  in  references  A. 10  and  A.  11. 
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The  advantages  of  the  FFP  are  that  it  includes  the  "nearfield",  and  it  can 
include  shear  propagation  effects  with  no  restrictions  on  the  shear 
velocity  relative  to  the  water's  congressional -wave  velocity.  Its 
disadvantages  are  that  the  procedure  is  not  easily  automated,  only  one 
source/receiver  depth  configuration  can  be  done  for  each  numerical  solution 
(see  Eqs.  A. 6,  A.  13  and  A.  14)  and  the  sampling  of  the  field  in  range  is 
restricted  by  Eq.  A. 12. 


A.  2  Normal-Mode  Solution 

The  alternative  to  a  direct  numerical  integration  of  Eq.  A. 6  is  to  expand  u 
into  a  complete  set  of  normal  modes: 

u(nx,ny;z)  =2  an(nx,ny)un(z),  (Eq.  a.15) 

where  the  un's  are  the  solutions  to  the  eigenvalue  equation 
d2u  (z) 

— ap—  +  lk2U)  -  knlun<z)  =  0  (Eq.  A16) 


that  satisfies  the  above-mentioned  boundary  conditions.  In  addition,  we 
require  the  u  (z)  be  bounded  as  z  -»  ».  The  normal  modes  u  (z)  form  a 
complete  orthonormal  set  that  satisfies  the  relation  ' 

1  p(z)un(z)um(z)dz  =  6nm,  (Eq.  A. 17) 

Jo 


where  the  density  p(z)  takes  its  appropriate  value  in  each  layer  and  6nm  is 

the  Kronecker-del ta  symbol.  The  spectrum  of  e:genvalues  consists  of  a 
discrete  part  and  a  continuous  part,  the  discrete  eigenvalues  occurring  in 
the  interval 


u)/c2  <  kn  <  max[u>/c(z)] ,  (Eq.  A.  18) 

where  c^  is  the  highest  speed  of  the  system.  In  the  present  treatment  we 

consider  only  the  discrete  eigenvalues,  since  in  general  the  continuous 
spectrum  makes  a  negligible  contribution  beyond  the  nearfield  of  the  source 
(and  requires  an  FFP-type  calculation  in  any  event). 


We  now  substitute  Eq.  A.15  into  Eq.  A. 6,  multiply  the  resulting  equation  by 
p(z)um(z),  and  integrate  over  z  from  0  to  giving: 


1  P(zo)un(zo) 
2n  q*  -  k2 


(Eq.  A. 19) 
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Substituting  Eqs.  A.  15  and  A.  19  back  into  Eq.  A.  5  we  obtain  an  integral 
representation  of  the  velocity  potential, 


<t>(x,y  ,z)  = 


p(zo}  r*  r  H  vji 

(W  J  -oo  n*  )-«  ny  \ 


W 

n2  *  kn 


x  exp[i(nxx  +  nyy)]  .  (Eq.  A. 20) 

We  evaluate  the  above  integral  by  choosing  a  path  about  the  poles  so  as  to 
lead  to  an  outgoing  wave  from  the  source  point  r  =  0.  Each  integral  in 
Eq.  A. 20  is  proportional  to  the  two-dimensional  plane-wave  representation 
of  the  zero-order  Hankel  function  of  the  first  kind  <A.12>,  and  therefore 
<Kx>y.z)  can  be  expressed  as: 


i))(r,z)  =  1  p(zQ)  Yj 
n 


u„(z.>,1(z)H<1)(k  r) 


(Eq.  A. 21) 


The  asymptotic  form  of  the  Hankel  function  can  then  be  used  to  obtain 


<t>(r,z) 


ip<zo>  -i„/4 
C  6 
(8/tr)^ 


E 


un<Zo>Un^ 


i  k  r 
n 


(Eq.  A. 22) 


In  addition  to  the  decay  of  the  field  due  to  cylindrical  spreading,  other 
loss  mechanisms  such  as  volume  attenuation  in  the  water  column  and  bottom 
are  included  in  Eq.  A. 22  because  the  eigenvalues,  k  ,  have  positive 
imaginary  parts  <A.13>,  thereby  resulting  in  an  exponential  attenuation  of 
each  normal-mode  term.  Equation  A. 22  gives  us  the  important  result  that 
the  field  at  a  depth  z  is  proportional  to  a  sum  of  the  products  of  normal 
modes  evaluated  at  the  source  and  the  receiver  depth.  The  normal  modes  are 
the  "natural  vibrations"  of  the  system  and  if  a  point  source  is  located  at 
the  null  of  a  particular  normal  mode,  that  mode  will  not  be  excited. 
Similarly,  if  a  point  receiver  is  placed  at  the  null  of  a  particular  mode, 
that  mode  contribution  to  the  total  field  will  not  be  sensed. 


In  analogy  to  the  FFP  procedure,  the  main  numerical  effort  for  the 
normal-mode  procedure  is  the  solution  of  the  eigenvalue  problem  defined  by 
Eq.  A.  16  and  the  boundary  conditions.  There  are  many  techniques  to  solve 
this  equation  <A.14,  A.15>  but  they  are  mainly  applicable  to  low-frequency 
or  shallow-water  propagation,  where  there  is  only  a  small  number  of  modes 
<A.16>.  However,  there  are  also  techniques  to  handle  deep-water 
high-frequency  propagation  using  normal  modes  <A.17>. 

The  advantages  of  the  normal-mode  procedure  are,  first,  that  once  Eq.  A.  16 
is  solved  we  have  the  solution  for  all  source/recei ver  configurations,  and, 
second,  that  the  whole  solution  procedure  can  be  highly  automated.  In 
addition,  the  normal-mode  procedure  can  be  easily  extended  to  slightly 
range-dependent  environments  using  the  adiabatic  approximation  where  mode 
coupling  is  neglected  <A. 18  to  A.23>;  numerical  methods  for  including 
mode-coupling  effects  are  present  areas  of  research  <A.24  to  A.27>.  The 
disadvantages  of  the  normal -mode  solution  are  that  conventional  procedures 
do  not  include  nearfield  effects  (the  exception  is  Stickler's  work  <A.28>, 
but  even  there  the  nearfield  is  evaluated  with  a  procedure  similar  to  the 
FFP  approach)  and  there  are  restrictions  on  how  one  can  treat  shear 
propagation  in  the  bottom. 
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Both  FFP  and  normal  modes  are  solutions  of  Eqs.  A.  5  and  A.  6.  The 
difference  between  the  two  is  that  the  normal -mode  method  restricts  the 
integration  to  horizontal  wavenumber  (q  is  a  continuously  varying 
horizontal  wavenumber,  whereas  k  are  the  discrete  horizontal  wavenumbers) 
in  the  interval  corresponding  Po  the  discrete  portion  of  the  spectrum 
defined  by  Ineq.  A. 18.  From  Eq.  A. 20,  we  see  that  the  integrand  has  poles 
at  n  =  k  .  Hence,  the  function  u  in  Eq.  A.  10  and  Xm  of  Eqs.  A.  13  and  A.  14 
should  have  poles  at  the  same  locations. 

Figure  A. 2  shows  a  plot  of  the  FFP  integrand  for  case  2A.  The  horizontal 
axis  is  the  wavenumber  (q)  and  the  vertical  axis  is  a  normalized  absolute 
value  of  the  integrand  amplitude.  The  vertical  dashed  line  separates  the 
discrete  spectrum  (defined  by  Ineq.  A. 18)  and  the  continuous  spectrum.  A 
blowup  of  the  discrete  spectrum  is  shown  in  Fig.  A. 3.  In  this  particular 
case  there  are  eleven  modes.  The  plot  shown  in  Fig.  A. 3  is  an  FFP 
integrand  and  the  asterisks  on  the  plot  are  the  locations  and  the 
amplitudes  of  the  modes  from  a  NM  calculation.  We  see  that  the  eigenvalues 
as  calculated  from  the  NM  program  (Table  A.l)  coincide  precisely  with  the 
peaks  (poles),  and  also  that  the  NM  amplitudes  correspond  to  the  amplitudes 
of  the  peaks.  (It  turns  out  that  this  one-to-one  correspondence  in  the 
amplitudes  is  because  there  is  virtually  no  loss  in  this  problem. 
Otherwise,  the  correspondence  would  not  be  as  precise  because  loss  shows  up 
in  the  FFP  calculation  as  widths  in  the  peaks;  nevertheless,  for  realistic 
losses  the  location  of  the  poles  would  be  the  same. ) 


TABLE  A.l 


MODAL  EIGENVALUES  FROM  SNAP 


Mode  no. 

IMiiliMB 

1 

0.10423 

2 

0.10386 

3 

0.10356 

4 

0.10328 

5 

0.10297 

6 

0.10261 

7 

0.10221 

8 

0.10182 

9 

0. 10142 

10 

0.10102 

11 

0.10064 

Figure  A. 4  displays  the  amplitudes  of  the  eleven  modes  plotted  as  a 
function  of  depth.  The  dashed  line  indicates  the  source/receiver  position. 
Notice  the  high  excitation  of  the  second  mode  and  the  low  excitation  of  the 
third,  sixth,  and  ninth  modes,  and  notice  the  one-to-one  correspondence 
with  the  FFP  integrand  shown  in  Fig.  A. 3. 

This  particular  example  shows  the  consistency  between  FFP  and  NM 
calculations.  In  the  next  section  we  show  the  results  on  transmission  loss 
when  discrete  and  discrete-plus-continuous  spectra  are  included  and  we 
compare  them  with  PE  calculations. 
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A. 4  Parabolic  Equation 

If  the  environment  varies  both  in  range  and  depth,  the  wave  equation  cannot 
be  separated  and  therefore  direct  numerical  integration  is  required.  At 
present  there  are  no  practical  methods  to  perform  this  direct  integration 
of  the  three-dimensional  wave  equation,  which  is  a  boundary- value  problem. 
An  alternative  approach  is  to  derive  an  approximate  wave  equation  that 
lends  itself  to  practical  numerical  solution.  We  now  outline  the 
derivation  of  such  an  approximation,  the  Parabolic  Equation  (PE)  <A.29>. 

The  velocity  potential  is  decomposed  as  follows: 

♦  =  <Kr,  2)  •  S(r) ,  (Eq.  A.23) 

and  we  substitute  <j>  into  Eq.  A.l  in  a  source-free  region: 


>l>  • 


i  as" 

r  8r 


+ 


S 


A.  24) 


That  is,  we  will  eventually  end  up  with  an  equation  that  allows  a 
"marching- in- range"  solution  and  we  will  have  to  initialize  the  solution  in 
some  way  (see  below).  We  use  the  notation  that 

k2  =  k2n2  ,  (Eq.  A. 25) 

where  n  is  an  "index  of  refraction"  equal  to  c  /c,  where  c  is  a  reference 
speed. 


Equation  A. 24  may  be  separated  into  two  differential  equations  by  setting 
the  terms  in  the  first  bracket  equal  to  -Sk2  and  the  terms  in  the  second 
bracket  equal  to  »|>k2  where  k2  is  the  separation  constant.  The  functions 
S(r)  and  tj»(r,z)  then  nave  to  satisfy  the  following  two  equations: 


a2s 

dr7 


+ 


A  +  1(2$  - 

r  3r  V 


(Eq.  A. 26) 


*  fr  +  kon2*  "  ko^  =  0> 


(Eq.  A. 27) 


The  solution  of  Eq.  A. 26  is  the  zero*'  order  Hankel  function,  ^(k  r), 
whose  asymptotic  form  has  been  given  in  Sect.  A.l.  Substituting  °the 
asymptotic  form  of  the  Hankel  function  into  Eq.  A. 27  and  making  the 
paraxial  approximation 


«  2k 


8ijt 

o  dr  ’ 


(Eq.  A. 28) 
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we  obtain 


(Eq.  A. 29) 


which  is  the  parabolic  wave  equation. 

The  paraxial  approximation  is  a  narrow-angle  approximation  It  implies 
that  the  rapid  range  dependence  of  Eq.  A. 23  is  included  in  S(r),  while  t|»  is 
a  function  varying  more  slowly  in  r.  An  approximation  to  solving  Eq.  A. 29 
is  to  assume  that  n  is  not  a  function  of  the  spatial  variables  but  is  a 
constant.  It  is  shown  elsewhere  <A.29,  A.30>  that  the  error  introduced  can 
be  made  arbitrarily  small  by  using  numerical  methods.  With  n  a  constant, 
we  can  Fourier  transform  41  with  respect  to  z, 

Hr, *)  =  ^  I  4>(r,z)e  iszdz,  (Eq.  A. 30) 

J  ~  00 

which  together  with  Eq.  A. 29  gives 

-s2i|>  +  2i k  +  k2(n2  -  1)  =  0.  (Eq.  A.  31) 

0  dr  0 


Equation  A.  31  is  a  first-order  differential  equation  with  constant 
coefficients  and  has  the  solution 


4>(r,s)  =  Hro,s) 


k2(n2  -  1)  -  s2 

. 


(r  -  ro) 


(Eq.  A. 32) 


where  the  initial  condition  at  r  must  be  specified.  The  field  as  a 
function  of  depth  is  the  inverse  transform  of  Eq.  A. 30 


4»(r,z) 


(n2  -  l)Ar 

s)  •  e 


e 


1&L  2 
2ko 


e1szds  (Eq.  A.33) 


where  Ar  =  r  -  r  . 

o 

By  introducing  the  symbol  for  the  fourier  transform  from  the  z-domain 
and  as  the  inverse  transform,  Eq.  A.33  may  be  written  as 

ik 


<Ji(r  +  Ar.z)  =  e 


-^(n2  -  l)Ar 


•/ 


’-1 


l^Ic2 

2ko 


J^Hr.z)^ 


(Eq.  A.  34) 
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Equation  A. 34  is  the  so-called  "split-step"  marching  solution  of  the 
parabolic  equation.  The  fourier  transforms  are  performed  using  an  FFT.  It 
is  the  solution  for  n  constant,  but  the  error  introduced  when  n  (profile 
and  bathymetry)  varies  with  range  and  depth  can  be  made  arbitrarily  small 
by  increasing  the  transform  size  and  decreasing  the  range-step  size  <A.29, 
A.  30>. 

The  parabolic  equation  is  not  a  boundary- value  equation  as  we  have 
numerically  formulated  it  above.  We  can  include  the  surface  boundary 
condition  by  taking  an  anti-symmetric  FFT  about  the  sea  surface  (z  =  0).  In 
practice  this  is  performed  by  taking  sine  transforms.  The  boundary 
conditions  in  the  bottom  are  simulated  by  including  the  discontinuity  in 
velocity  in  the  sound-speed  profile.  There  are  methods  to  also  include  the 
density  discontinuity  <A.29>.  The  radiation  condition  as  z  goes  to 
infinity  is  simulated  by  requiring  the  field  to  exponentially  tail  off  for 
large  values  of  z  beyond  which  there  would  not  be  any  significant  acoustic 
interaction. 

As  mentioned  above,  the  PE  method  requires  an  initial  starting  solution. 
Two  methods  have  been  used  for  describing  a  point  source.  The  first  method 
is  to  initialize  the  field  with  a  set  on  normal  modes  descriptive  of  the 
point  source  in  the  starting  environment.  This  would  not  include  the 
continuous  portion  of  the  spectrum  (see  Sect.  A. 3),  but  for  long-distance 
propagation  this  approximation  is  adequate.  A  second  approach  has  proved 
to  be  simpler  and  as  effective.  The  point  source  is  approximated  by  two 
gaussians  that  are  anti -symmetric  about  the  sea  surface,  thereby  auto¬ 
matically  including  the  pressure- release  boundary  condition  at  the  surface. 
It  is  this  last  technique  that  has  been  used  in  this  report  and  we  have 
found  that  not  only  is  this  a  good  approximation,  but  we  have  also  shown 
that  by  using  this  technique  part  of  the  continuous  spectrum  is  included. 

In  Sect.  A. 3  we  compared  FFP  and  NM,  using  test  case  2A  as  an  example.  We 
now  look  at  transmission  loss  from  the  point  of  view  of  discrete,  discrete- 
plus-continuous,  and  PE,  which  does  not  obviously  distinguish  between 
regions  in  the  spectrum.  Figure  A.  5  shows  test  case  2A  plotted  from  0 
to  6  km.  Note  that  the  PE  tracks  the  FFP  results  in  the  nearfield 
indicating  that  at  least  part  of  the  continuous  portion  of  the  spectrum  is 
included  in  the  PE  calculation.  Figure  A. 6  for  the  full  range  of  20  km, 
shows  how  the  three  model  results  converge  in  the  farfield;  recall  that  the 
NM  calculation  does  not  include  the  full  nearfield  contribution.  This 
particular  example  clearly  illustrates  the  consistency  and  inter¬ 
relationship  between  the  three  models. 

The  advantages  of  the  PE  are  that  it  handles  a  range-dependent  environment 
and  gives  the  acoustic  field  in  the  entire  water  column  without  additional 
computational  effort.  Its  disadvantages  are  that  the  procedure  is  not 
easily  automated,  and  it  is  practical  only  for  low-frequency  propagation 
since  computation  times  increase  with  frequency  squared.  Moreover,  there 
is  no  straightforward  way  of  handling  shear  propagation  in  the  bottom. 


A. 5  Ray  theory 

This  report  is  concerned  mainly  with  wave  theory;  nevertheless,  for 
completeness,  we  include  a  brief  description  of  ray  theory.  In  this  case 
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we  assume  a  solution  of  Eq.  A.l  (with  right-hand  side  equal  to  zero)  as 

4>  =  >|i(x,y,z)  *  eiS(x’y’z).  (Eq.  A.35) 

S(x,y,z)  is  a  phase  function  that  includes  rapid  variations  as  a  function 
of  range,  and  <J»(x,y,z)  is  a  more  slowly  varying  envelope  function  in  which 
geometrical  spreading  and  loss  mechanisms  are  included  (in  the  PE,  S 
contains  the  cylindrical  spreading  factor).  Substituting  Eq.  A.35  into  the 
wave  equation  and  separating  real  and  imaginary  parts,  we  obtain 

|  V2i|i  -  (VS)2  +  k2  =0,  2(V<j>  *  VS)  +  >|)V2S  =  0.  (Eq.  A.  36) 


We  now  make  the  geometrical -acoustics  approximation 

V2i|i  «  k2,  (Eq.  A.  37) 

that  is,  the  amplitude  of  the  phase  function  varies  slowly  in  range  with 
respect  to  wavelength.  Substituting  Eq.  A. 37  in  Eq.  A. 36  gives  the  eikonal 
equation, 

(VS)2  =  k2.  (Eq.  A. 38) 


The  trajectory  of  the  rays  is  perpendicular  to  the  surfaces  of  constant 
phase  (wavefronts),  S,  and  is  expressed  by 


(Eq.  A. 39) 


where  £  is  the  arc  length  along  a  ray  and  X  is  the  coordinate.  It  can  be 
shown  that  the  direction  of  the  average  energy  flow  is  along  these 
trajectories  and  the  amplitude  of  the  field  at  any  point  can  be  obtained 
from  the  density  of  these  rays;  formally,  having  solved  for  S,  the 
amplitude  is  obtained  from  solving  the  second  part  of  Eq.  A. 36.  We  also 
mention  here  that  corrected  ray  theory  assumes  that  i|»  is  a  function  of 
frequency  and  an  expansion  in  powers  of  inverse  frequency  is  made,  the 
leading  term  being  the  infinite-frequency  solution  with  the  additional 
terms  being  corrections  from  the  infinite- frequency  solution. 

The  advantages  of  ray  theory  methods  are  that  the  computations  are  rapid 
and  that  ray  traces  give  a  very  physical  picture  of  the  acoustic  paths. 
The  disadvantage  is  that  ray- theory  is  an  infinite- frequency  approximation 
and  therefore  does  not  include  diffraction  and  other  wave  effects.  This 
shortcoming  also  prevents  ray  theory  from  adequately  describing  significant 
bottom  interaction  and  low-frequency  ducted  propagation. 
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FIG.  A.  2  FFP  INTEGRAND  PLOT  FOR  CASE  2A  (full  spectrum). 
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FIG.  A. 3  FFP  INTEGRAND  PLOT  FOR  CASE  2A  (discrete  spectrum). 
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FIG.  A.  4  MODE -AMPLITUDE  FUNCTIONS  FOR  CASE  2 A 
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FIG.  A.  5  INTER-MODEL  COMPARISON  FOR  CASE  2A  (near  field) 
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FIG.  A. 6  INTER-MODEL  COMPARISON  FOR  CASE  2A  (far  field) 
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APPENDIX  B 

DESCRIPTION  OF  THE  FOUR  COMPUTER  MODELS 


B.l  The  FFP  Model 

The  FFP  model  <B.l,  B.2>  is  a  complete  numerical  solution  of  the  range- 
independent  wave  equation  and  hence  includes  the  continuous  part  of  the 
spectrum  and  the  nearfield  (distances  greater  than  a  wavelength).  The 
model  used  in  this  study  was  developed  at  Columbia  University  <B.2>.  It 
includes  propagation  of  both  compressional  and  shear  waves  and  is  therefore 
also  suitable  for  seismic  studies.  In  the  original  version,  the  depth- 
dependent  part  of  the  wave  equation  was  solved  using  the  Thompson-Haskell 
matrix  method,  where  the  water  column  is  divided  into  isovelocity  layers. 
An  updated  version  <B.3>  allows  constant  gradient  in  k2  (equivalent  to 
1/c2  constant  gradient)  and  therefore  employs  Airy  functions  in  the 
Thompson-Haskell  matrices  instead  of  trigonometric  functions.  The  inputs 
are  the  same  as  those  of  the  normal-mode  model  (SNAP).  There  are  two 
possible  outputs.  One  output  is  the  "integrand"  plot,  which  is  essentially 
a  plot  of  energy  versus  wavenumber  and  therefore  has  maxima  corresponding 
to  the  normal-mode  eigenvalues  in  the  discrete  part  of  the  spectrum;  the 
maxima  in  the  continuous  portion  of  the  spectrum  correspond  to  the 
so-called  "virtual  modes".  The  second  output  is  loss  versus  range.  It 
should  be  noted,  however,  that  the  present  program  requires  a  complete  new 
run  for  a  change  of  source  or  receiver  depth;  hence,  loss  contours  as 
functions  of  range  and  depth  are  not  practical,  as  compared  with 
normal-mode  or  PE  methods. 


B.  2  The  SNAP  Model 

The  SNAP  Model  <B.4>  is  a  normal-mode  model  based  on  a  program  originally 
developed  at  the  US  Naval  Research  Laboratory  <B.5,  B.6>;  this  program 
solves  the  eigenvalue  problem  by  direct  numerical  integration  of  the  depth 
dependent  equation.  Computation  time  for  some  of  the  key  subroutines  has 
been  reduced  and  the  program  has  been  restructured  to  run  interactively  on 
a  UNIVAC  EXEC-8  system.  The  model  allows  for  slight  range  dependence  by 
employing  the  adiabatic  approximation.  SNAP  was  originally  designed  for  a 
shallow-water  environment  but  there  is  now  a  version  that  handles 
high-frequency  deep-water  situations. 

Environmental  inputs  are: 

•  arbitrary  sound  speed  profile  as  function  of  depth  (multiple 
profiles  for  range-dependent  adiabatic  computations)  in  the  water 
column, 

•  density,  attenuation,  and  compressional-velocity  profile  of  sediment 
layer, 

•  density,  shear  velocity,  compressional  velocity,  shear  attenuation, 
and  compressional  attenuation  of  the  basement. 
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At  present,  the  output  options  are: 

•  loss  vs  range, 

•  loss  vs  depth, 

•  depth- averaged  loss  vs  range, 

•  depth-averaged  loss  vs  frequency, 

•  contoured  loss  vs  frequency  and  range, 

•  contoured  depth-averaged  loss  vs  frequency  and  range, 

•  modal  group  velocity  vs  frequency, 

•  modal  phase  velocity  vs  frequency, 

•  mode  function  vs  depth, 

•  phase  of  field  vs  depth, 

•  intensity  of  field  vs  arrival  angle, 

•  sound  speed  vs  depth. 

The  SNAP  model  has  been  subjected  to  extensive  inter-model  comparisons 
<B.7>  as  well  as  a  thorough  testing  against  experimental  data  <B.8,  B.9>. 


B.  3  The  PAREQ  Model 

PAREQ  <B.10,  B.ll>  is  a  parabolic  equation  model;  the  computer  program  used 
in  this  study  is  a  modified  version  of  the  one  developed  by  the  Acoustic 
Environmental  Support  Detachment  (AESD)  of  the  Office  of  Naval  Research 
<B.12>.  This  model  not  only  handles  a  variable  sound-speed  profile  in 
depth  and  range  but  also  allows  the  bottom  depth  and  bottom  structure  to 
vary  in  range.  The  bottom  is  characterized  by  a  compressional -velocity 
profile,  density  (the  density  discontinuity  is  smoothed  using  a  hyberbolic 
tangent  function  as  suggested  in  <B.10>)  and  attenuation,  which  is  included 
by  using  a  complex  sound  speed.  The  second  layer  of  the  bottom  has 
constant  acoustic  properties.  There  are  two  options  for  the  initial  field: 
gaussian  source  or  normal  modes.  The  PAREQ  model  is  not  only  resident  on  a 
UNIVAC  1100/60  but  also  runs  on  an  HP  21MX  computer.  The  model  has 
recently  been  used  to  study  mode  conversion  and  mode  cut-off  phenomena  in 
range-dependent  environments  <B.13>. 


B.  4  The  GRASS  Model 

GRASS  <B.14>  is  a  range-dependent  ray-trace  program  originally  developed  at 
the  US  Naval  Research  laboratory  for  a  CDC  3800  computer.  This  version 
runs  on  a  UNIVAC  system;  the  only  significant  difference  between  the  model 
used  in  this  study  and  the  original  is  that  the  input  procedure  has  been 
simplified. 
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